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Abstract 

A method to sum over logarithmic potential in 2D and Coulomb potential in 3D with periodic 
boundary conditions in all directions is given. We consider the most general form of unit cells, the 
rhombic cell in 2D and the triclinic cell in 3D. For the 3D case, this paper presents a generalization 
of Sperb's work [R. Sperb, Mol. Simulation, 22, 199-212(1999)]. The expressions derived in this 
work converge extremely fast in all region of the simulation cell. We also obtain results for slab 
geometry. Furthermore, self-energies for both 2D as well as 3D cases are derived. Our general 
formulas can be employed to obtain Madelung constants for periodic structures. 
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I. INTRODUCTION 



It has become a common practice to employ numerical simulations in the study of phys- 
ical problems, which are difficult to solve analytically. Since it is not possible to simulate 
realistic physical systems, containing ions of the order of Avogadro number, one usually 
works with a very small system. For small systems, containing a few hundred to a few 
thousand charges, boundary effects become relatively pronounced, especially if the nature 
of interaction is long range. To avoid this problem, periodic boundary conditions (PBC) are 
usually employed. In many simulations, the nature of interaction is such that the potential 
satisfies the Poisson equation. For example, a logarithmic interaction in two dimensions 
(2D) and a Coulomb potential in three dimensions (3D) both satisfy the Poisson equation in 
2D and 3D respectively. We refer to a potential which goes as r^' d+2 ^ in a d > 2 dimensional 
isotropic space as a Coulomb type potential. The Coulomb type potentials fall under the 
category of long range potentials. In fact, in a d > 2 dimensional space, any interaction 
which goes as r~ a , where a < (d — 1) is known as a long range interaction. The reason 
being, while the potential decays as r~ a , the volume element goes as r^" 1 ). As a result, in a 
periodic system, even charges located at infinity give rise to a finite contribution to energy 
and forces, which cannot be neglected. We consider Coulomb type of potentials in 2D and 
3D in this paper. 

To derive a formula for interaction between two particles with PBC imposed, one has to 
consider the interaction of a particle with periodic repetitions of itself, as well as that of the 
second particle. Interaction energy of a particle with its own periodic repetitions is termed 
as the self-energy. Determination of self-energy is important in simulations where the size 
of simulation box may change during the simulation. For example, such a case arises in an 
isobaric Monte Carlo simulation. The aim of this paper is to consider the kind of interactions 
mentioned above for the most general type of unit cells in 2D and 3D. We consider a rhombic 
unit cell in 2D and a triclinic cell in 3D, with origin lying at the bottom left corner of the 
unit cell. The unit cell contains a number of ions, which interact via Coulomb type potential, 
satisfying the Poisson equation in their respective dimensions. The unit cell repeats itself 
in all directions under PBC. Hence, the interaction of a particle located at r with another 
particle located at the origin includes, apart from the direct interaction between the two 
particles, the interaction of the first particle with all periodic images of the second particle. 
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These periodic images of the second particle are located at lattice vector sites given by 
1 = ma + nb + pc, where m, n and p range from — oo to +00. Also, the particle interacts 
with its own images located at r + 1, where 1 is defined as above. Thus, if we have N charges 
qi in a charge neutral unit cell, then the Coulomb energy may be written as 



where a prime indicates that n = term is to be excluded for the case when % = j. 
The series in Eq. fjl.lj) is a conditional series. This series can be summed up to any value 
depending on the order in which the terms of the series are grouped. Therefore, a summation 
convention has to be specified based on the physical nature of the problem in mind. 

The conditional series mentioned above may be evaluated by introducing background 
charges in a way that the total background charge adds up to zero. Imposing background 
charges in this way leads to well defined ways of summing the conditional series. However, 
results of the summation of conditional series may still differ in view of the method employed 
to impose background charges, as the background charges may have a structure of their own. 
For example, background charges in a 3D system with PBC may be imposed in the following 
two ways. A charge q and all its periodic repetitions under the PBC may be viewed as a set 
of layers along an axis of the unit cell. In order to impose background charge on this system, 
we may assume that all these different layers are charge neutral separately. Thus, we may 
assume that for a layer composed of charge q and its periodic images, one has an additional 
uniform charge density of —q/a, where a is the area of the 2D unit cell. Thus the overall 
charge contained in each 2D cell of the layer is zero. Another charge q' present in the system 
will interact with the set of charges q as well as the neutralizing background surface charge 
with charge density. However, it can be shown that introducing these uniform background 
charge sheets leads to some unwanted terms, as the sheets have a structure of their own. 

However, there is a better way of imposing background charges, without introducing 
any structure of the background charges themselves. This can be achieved by distributing 
a uniform 3D charge on the grid made out of charge q and its images. The neutralizing 
background charge now has a uniform charge density of —q/V, where V is the volume of 
the unit cell. This volume charge adds up to zero at any point due to the overall charge 
neutrality condition and thus does not introduce any artificial structure, such as the uniform 
sheets in the previous case. 




N 



(1.1) 
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The results of the two prescriptions suggested above differ by two terms. The first term 
depends on the square of the component of the dipole moment, along the direction of 
layering, of the original charges contained in the unit cell. The second term depends linearly 
on the distance between the pair of charges along the direction of layering. 

In this paper, we adopt the second procedure. Using the results derived here, it will be 
easy to establish connection between the results of two summation conventions mentioned 
above. Introduction of neutralizing charge background in the form of a uniform cloud leads to 
only the intrinsic parti of potential energy and this technique has previously been employed 
by LeknerJ» and Sperb 2 . It is important to know that the two procedures mentioned above 
still do not lead to the correct energy of a collection of charges interacting under the PBC, 
if one wants a limit of spherical means 2 -. De Leeuw et. al 3 have shown that for 3D case, an 
extra term depending on the total dipole moment of the unit cell has to be added to get the 
correct energy of charges. For the 2D case the correction term turns out to be zero. 

With the help of discussion above, the energy of N particles contained in a unit cell 
with periodic boundaries and interacting through a Coulomb type potential in 3D can be 
expressed as, 

Here the charges are denoted by q^s and their positions in the unit cell by r-j's and 1 < i < N. 
The last term in Eq. (jl.2|) is the dipole term introduced by De Leeuw et. al 3 . For the 2D 
case one has only the first two terms on the right hand side. Our aim in this paper is to 
obtain expressions for G(r) and G 8e n in 2D and 3D. 

Before proceeding further, we briefly discuss three main approaches in use to obtain 
Coulomb interaction with periodic boundaries. These three approaches are due to Ewald^, 
Lekneri^ and SperfA The Ewald method was developed eighty years ago in connection with 
the evaluation of Madelung constants. This method, in spite of its shortcomings, is still very 
much in use. The method proceeds by breaking the original summation in two parts. One 
of these sums is carried out in real space and another one in Fourier space. This splitting of 
summation depends on a real parameter which has to be chosen judiciously, failing which 
the series in real and Fourier space might converge very slowly. In general, to calculate a 
pair-wise interaction it usually requires a few hundred terms involving complementary error 
functions. 
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An alternative to the Ewald method was given by Lekner—. This method involves an 
evaluation of a few dozen terms if the position vector r is not very small. If r tends to 
zero this method converges slowly. The problem of convergence was fixed later by Lekner 
in another paper—, following Sperb's work 6 . Lekner method has not been generalized to a 
triclinic cell yet. Though, it is possible to generalize Lekner's work for a triclinic cell, here, 

to obtain results for a triclinic cell. 

Among the latest advances on Coulomb sums is by Gr0nbech- Jensen^ in 2D and Sperfj 2 - 
in 3D. Sperb's results are similar to that of Harris et a/.— and Crandall et a/.—. A major 
advantage of Sperb's work is that it can be employed to get N\n(N) scaling in timeii, where 
N is the number of ions present in the system. On the other hand, with Ewald summation 
method, one can get only [N ln(iV)] 3//2 scaling 12 . In this paper, our aim is to generalize 
Sperb's work 2 to a triclinic cell. The method given in this work will contain Sperb's result 
in a simplified form as a special case. Also for the first time, an alternative to Ewald's 
technique will be given, which can be applied to the most general kind of unit cell in a 
computer simulation, a triclinic unit cell. We will also discuss scaling of Nln(N) that may 
be achieved with the use of formulas developed here. 




FIG. 1: A schematic diagram explaining the set of for a 2D rhombic cell. A charge located at B 
interacts with another charge located at the origin A as well as its periodic images located at A's. 
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II. LOGARITHMIC INTERACTION IN 2D 



An excellent method to sum over Coulomb type potential (logarithmic interaction) in 2D 
was given by Gr0nbech- Jensen^. Another alternative was provided by Tyagi et alM in a 
recent paper. The problem with the later approach is that the lattice sum does not converge 
when the two charges are close together within the unit cell. This problem will be addressed 
here and formulas will be modified in such a way that the convergence is achieved for even 
those cases where charges are close to each other. Thus we will obtain a result which is 
different from Gr0nbech- Jensen but still as efficient. 

We consider a rhombic cell with periodic boundaries along the x and y directions. A sketch 
of the cell is shown in Fig. 1. A particle, located at position r, interacts logarithmically 
with charges located at the vertices of a rhombic grid. A formula was developed in Refill] to 
compute this sum. We sketch a portion of that derivation here for the sake of completeness. 
Consider the Poisson equation in 2D: 

V 2 G(r) = -27T 8(r + 1) + (2.1) 

i 

The second term on the right hand side amounts to the presence of a neutralizing background 
charge. The solution of Eq. (|2.1|) is given by 

n( \ 2n v ex P(*Q " r ) 1 \ (oo\ 

G(r) = Tj^e h -cFHT -?)• (22) 

where l\ and l 2 denote the lengths of the sides of the rhombic cell and 

r = r 1 e 1 + r 2 e 2 , Q = nibi + n 2 h 2 . (2.3) 

Here, < r\ < l\, < r 2 < h arid ei and e 2 represent the unit vectors along the axis of the 
rhombic cell, with ei.e2 = cos#. We have also introduced an infinitesimal parameter £. The 
sum over Q runs over all reciprocal lattice vectors spanned by 

2tt 

h i = > ■ 2n ( e i ~ ejCOSfl), (2.4) 

Li sin 6 

for = (1, 2), (2, 1) and nx and n 2 are integers. Introduction of an infinitesimal parameter 
£ as in Eq. (J2.2j) implies assumption of the presence of a neutralizing background charge. 
Thus, here a charge q located at (x, y) in the unit rhombic cell interacts with charges q' 
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located at the origin and all other vertices of the grid. The charge q also interacts with a 
uniform layer of background charge superimposed on the grid of q' charges such that the 
charge density is —q'/a, where a = Z 1 Z 2 sin^, is the area of the unit cell. From now onward 
we will always assume that a final limit £ — ► is to be taken. Using the value Q from 
Eq.Q and Eq.flHJ) in Eq.fl22J we obtain 



Q1T1 ft °° °° eX P 

gw-^E E 



sin# 1 



where < rj/ij < 1, u = and we have redefined infinitesimal parameter £ for sake of 
calculations. We now evaluate the sum 



exp ( ^7^2 



£ ^ ^ — ; 2c2 (2.5) 



IF 



H E 



exp I ij L n 2 r 2 



n 2 a 2 — 2n\n 2 a cos + ni + a 2 f 2 

n2=— 00 1 z 



,2 , ;0 _ , + . exp(-i27r/3 n J sinh [7^2] + sinh [27r7 m (l - t 2 ) 



nil exp(«27r/? ni t2) 



7 m [cosh(2vT7 ni ) - cos(27r/3 n J] 
where t 2 = r 2 /l 2 , 



(3 n = na cos 0, 7„ = a\/(n 2 sin 2 + £ 2 ), (2.6) 
and we have used the identity (here a < 2n) 



exp (ina) ir exp [if3 (a — 2n)] sinh (7a) + exp (i(3a) sinh [7 (2n — a)] 
nf^oo ( n ~ I 3 ? + 7 2 ~ 7 cosh (2tt7) - cos (2tt/3) 

which is derived in Appendix A. The sum defined in Eq. (J2.5|) can be written as 



a sin 9 + °° 

G(r) = £ ex P [*2vr (n*i + /3„* 2 )] (2.8) 

n=— 00 

exp(— 2nf3 n ) sinh(27T7 Tl t2) + sinh [27T7 n (l — t 2 )} sin 9 1 
7 n [cosh(27T7„) - cos(27r/5 n )] 2-no £ 2 ' 

where ti = T\/l\. Separating out the term corresponding to n = 0, we obtain 
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a sing / sinh(27rft 2 ) + sinh [2<(1 - t 2 )] \ sing 1 
GW " 2 V e [cosh(27rO - 1] J 2*a? { ' 

a sing exp(-27r/3 ra ) sinh(27r7 n t 2 ) + sinh [27T7 W (1 - t 2 )\ 

2 Z_^ ex Pi Z7m J 7n [cosh(27r 7n ) - cos(27r/5 n )] 

where £ = ti + t 2 a cos g and a prime on the summation sign indicates that the term corre- 
sponding to n = is not to be included. Taking the limit £ — > 0, one obtains 



V £cr [cosh(2<cr) - 1] na 2 £ 2 J 3 V 2; V ; 

Thus we obtain the following expression for G(r) 



. . cr sin g it , 9N , 

G(r) = - (1 - 6t 2 + 6t 2 2 ) (2.11) 

o- sin g >A exp (-i2?r/? n ) sinh (27T7 n0 t 2 ) + sinh [27T7 n0 (l - t 2 )] 

2 Z^P^ 7 ™; 7n0 [cosh(27r 7n0 ) - cos(2vr/?„)] 

where 7 n o = cr |nsing|. Due to the symmetrical nature of the unit cell, it suffices to look 
at only that part of the unit cell, which corresponds to < t\ < 0.5 and < t 2 < 0.5. 
Eq. (j2.11j) fails to converge fast enough when £ 2 — * 0. This problem can be easily fixed as 
follows. We add and subtract the following term from Eq. (|2.11|) 

h (t, h) = lfl eXP( ~ 27r7ra ° t f f eXp(2?rmt) . (2.12) 

The quantity h (t,t 2 ) can be easily evaluated by carrying out the sum in Eq. ()2.12j) analyt- 
ically. Using the identity 



g exp(-n|a|)cos(n6) = _ 1 ^ ^ ^ _ ^ (&)] _ lnj2) + Q (2 ^ 

one obtains 



n 2 2 2 

n=l 



h(t,t 2 ) = — -In (cosh[27rf 2 0"sing] — cos[27rt]) 
In (2) 27rcrt 2 sin g 



(2.14) 



S 



Thus we can write 



. . a sin 9 it , 9X , 

G(r)=^— - (1 - 6* 2 + 6*1) (2.15) 

1 , , , ,„ s , , 2nat 2 sm9 In (2) 

- - In (cosh [a sin (0) 2vrt 2 ] - cos[2?rt]) ' 



2 2 2 

1 >A f exp(-i27r/3 n ) sinh(27r7 n0 t 2 ) + sinh [27T7„ (1 - t 2 )) 
x - > exp (izirnt) < ; — — T — — : ; — 

2 ^ V ; 1 \n\ [cosh(27r 7n o) - cos(2vr/3 n )] 

exp(-27T7 n0 *2) 



\n\ 

After some effort, the above equation can be written as 

G(r) = (l + 6*1) - ]- In (cosh [2vra sin (0) t 2 ] -cos[2vrf]) (2.16) 

In (2) + °° 

h {cos [27r (n* — /3 n )) sinh(27T7 ra0 *2) + cos(27mt) 

n=l 

x [exp (-27T7 n0 t 2 ) cos (2vr/3 n ) - exp (-2vr7 n0 ) cosh (27T7 n0 i 2 )]} / 
{n [cosh(27T7 n0 ) - cos(27r/3 n )]} . 

Equation (|2.16|) converges extremely fast for all values of < t 2 < 0.5. However, to achieve 
better convergence the sides of the rhombic cell should be labelled such that a = l 2 jl\ > 1. 
Now, an expression for the self-energy can be easily obtained by taking the limits t\ — > 
and t 2 — ► and subtracting 

g (r) = -^ln (r\ + r 2 + 2rir 2 cos6 l ) , (2.17) 

one obtains 

G s 2 f = lua(G(r)-g(r)) (2.18) 
o sin 7r ^ f sinh(27T7 n o) 



2 3 \n [cosh(27T7 n0 ) — cos(27r/3 n )] n, 

- lim Q In [(2vrt 2 a sin 0) 2 + (27rf) 2 ] - i In [r 2 + r 2 + 2rxr 2 cos 0] ^ 

= - in (2./,) - v ( e ^(- 2 ^°) ~ cos( ^;M . 

2 3 v 17 ^ V^[ c osh(27T7 n0 ) -cos(2vr/3 n )]y 
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III. COULOMB INTERACTION IN 3D 



The Poisson equation to be solved in this case is 

■+- n + 

v 



Air 

V 2 G(r) = -47r^<5(r + l) + -, (3.1) 
i 



where 

V = hhh[e i .{e j x e fc )] (3.2) 

stands for the volume of the unit cell. The last term in Eq. (j3.1|) amounts to the presence of 
uniform background charge. The solution of Eq. ()3.1jl is given by 

ow-v&yZ-vTe—e)' (3 ' 3) 

where 

r = nei +r 2 e 2 + r 3 e 3 ,, Q = mbi + n 2 b 2 + n 3 b 3 , (3.4) 

where Q runs over all reciprocal lattice vectors spanned by 

2tt e,- x ej. , 

bi = -, * v (3.5) 

for all cyclic permutations of (i,j,k) and ni,ri2 and n 3 range from — oo to +oo. Using 
Eas. (f3~3]) . (j3~l|) and (jS3J) we obtain 



4tt ^ ex P i2n \ ni h +n2 h +n 3h) 

Vb\ n f^ ni ("J + n\c 2 2 + nf en + 2nin 2 c i2 + 2n 2 n 3 c 23 + 2n 3 n!C 3 i + £ 2 ) 
4vr 1 



where < rjli < 1, 

c ^ = trir l<i,i<3, (3.7) 

b 3 .b 3 

and we have, as before, introduced an infinitesimal parameter £ in the denominator and 
subtracted a counter term from the whole sum due to the presence of a uniform background 
charge. We now evaluate the sum 

1 ^ ex P fe^a) 
L(n\, n 2 r 3 , £) = — > — -5 5 5 — . 

ix ' rig + n 2 c 23 + nfci3 + 2nin 2 ci 2 + 2n 2 n 3 c 23 + 2n 3 nic 3 i + £ 2 ) 

723 = — OO 

(3.8) 
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This sum can be obtained easily and the result is 



where 



L (711,712,7*3,0 (3.9) 

Co a + \ exp H 271 "/^,^] sinh(27T7„ lin2 t 3 ) + sinh [2ir-f nun2 (1 - £ 3 )] 
= exp{t27rfJ nun2 t 3 ) = — — — — r= 

7ni,n 2 COSh(27T7n lin2 ) - COS(27T/^ niin2 ) 



U = y for i = 1,2 and 3, (3.10) 



Pn u n 2 = -7ilC 3 i - n 2 C 3 2, (3-H) 



and 



7m,na = [ n 2 c 22 + n'cn + 2n 1 n 2 c 12 - (711C31 + n 2 c 32 ) 2 + £ 2 ] 1/2 . (3.12) 
Plugging the value of L(ni,n 2 r 3 ,0 from Eq. (|3.9|) in Eq. ()3.6|) we obtain, 



An 1 

G ( r ) = 7/^2 ex P [i27r(raiti + n 2 t 2 )} L(m, n 2 r 3 , - T/52 72 ( 3 - 13 ) 

3 ni,n.2 ^ 
4-7T 2 ' 

exp [z27r(riiti + n 2 t 2 )] L(n 1 ,n 2 r 3 ,^) 



ni,n2 



f=0 



4ff 2 7T , 9N 



W|3 

where a prime over summation sign implies 77,1 and n 2 cannot both be zero simultaneously. 
We have also separated out the term corresponding to n\ = and n 2 = in Eq. (|3.13| ) and 
taken the limit £ — > 0, which results in cancellation of the diverging factor 47r/ (Vfe 2 ^ 2 ). Eq. 
()3.13j) is one of the main results of this paper. It is easy to see that the sum defined in 
Eq. (j3.13j) fails to converge fast enough as ti tend to zero. In fact, towards large values of 
7„ 1>n2 , the quantity L defined in Eq. (l3.13|) goes as exp (— 27rj nitn2 t 3 ) and if t 3 is small, this 
convergence may be very slow. As before for the 2D case, we only concentrate on that part 
of unit cell which corresponds to < £, < 0.5. To transform Eq. ()3.13|) in a form, which 
converges even for small values of ti, we need to separate out a term which corresponds to 
slab geometry. By slab geometry we mean a situation which is obtained by sending one of 
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the sides of the unit cell to infinity. We write 

+ G 2 (r) (3.14) 



G(r) = —2 ^2 ex P [«27r(ni*i + n 2 t 2 )} B (n x , n 2 r 3 , £) 



J ni,rt2 



5=0 



+ ^6|3 ( 1 + 6f ")' (3 ' 15) 

where G 2 corresponds to the slab geometry case and is given by, 



= ™| 2. — ( 3 - 16 ) 

47T 2 7T 

- TTfef 3 (6t3) ' 
and 5 is defined as 

o, c\ t( c\ exp(i27r/? niin2 t 3 )exp(-27r7 niin2 t 3 ) 

5(ni, n 2 r 3 , = L(m, n 2 r 3 , £) . (3.17) 

Tni,H2 

The result in Eq. (|3.16|) represents Coulomb interaction with open boundary condition along 
the r 3 direction and periodic boundaries along the t\ and r 2 directions. G 2 can be obtained 
from G by taking the limit Z 3 —>■ oo and dropping a constant term. We also note above that 
both (/3 niin2 t 3 ) and (7 ni ,n 2 ^3) are independent of l 3 , when £ — > 0. 
The term in Eq. ()3.17|) . can be written as 

B(m,n 2 r 3 ,0 (3.18) 
exp (z27r/3 ni)n2 t 3 ) [cosh [2ti (iP„ un2 - ln u n 2 h)] ~ ggH^wJ cosh(27r7 ni)n2 t 3 )] 

7m,n 2 [cOsh(27T7 niin2 ) - COs(27T/5 niin2 )] 

It can be easily seen that, for large 7 ni ,n 2 j the slowest decaying term on the right hand side of 
Eq. ()3.18| ) goes as exp [— 27T7 ni n2 (l — 1 3 )]. So, the fastest convergence now occurs for t 3 = 
and slowest for t 3 = 0.5. But even this 'slowest' convergence for t 3 = 0.5, amounts to an 
extremely fast exponential convergence of exp(— vr7 nijn2 ). 

Essentially now the whole problem has reduced to a fast evaluation of G 2 in Eq. (|3.16|) . 
We take up this case now. As G 2 (r) fails to converge fast enough for small separations, we 
break the sum in Eq. (j3.18|) into two parts 

E = E + E • ( 319 > 

"l,«2 ni,n 2 ni,n2=0 

12 



Note that in our notation Y' n an d Yin' represent the same thing. Thus G 2 (r) can be written 

as 

G 2 {t) = G' 2 (t)+G w {t), (3.20) 

where 

r>'(\ ^ ST co \ ( exp(2vrm 1 x 1 )exp(-27r7 niin2 t 3 )\ 

G 2i r ) = T7u 2^ exp ( t27Tn ^) 2^ z • ( 3 - 21 ) 



W>3 V 7ni,n 2 

J n2 \m=— oo ' 



and 



4tt 2 >A exp [i2Tm t (-c 31 t 3 + t{)] exp(-2vr7 niin2 t 3 ) 



4tt 2 



n 2 =0, £=0 



f 27rt 3 . 



(3.22) 



To further transform Eq. ()3.22|) . we express 7„ 1 ,n 2 as 

1 /2 

Tm.na = K ( c n _ 40 + n 2 ( c 22 - c 23 ) + 2n x n 2 {c X2 ~ ci 3 c 23 ) + £ 2 ] (3.23) 

r ~ A 2 ~ Nl/2 
n\8 + n 2 a + r?|& 2 



where we have put £ = and 

T I 2 \l/2 ~ (C12-C13C23) , . 

d=(cu-c 13 J , a= ~ , (3.24) 



1 /2 

£ = [(Cll - C 2 3 ) (C 22 - C^) - (C12 ~ Ci 3 C 23 ) 2 ] 

As the convergence of Eq. (j3.14j ) crucially depends on the value of 7ni,n.2> ^ is helpful at this 
point to note that the minimum value of 7 ni) „ 2 , when both n\ and n 2 are integers such that 
they cannot both be zero simultaneously, is given by 

7 2 lin = 6 2 min(l,-^-). (3.25) 
y a 2 + b 2 J 

We note that 7 min depends upon the geometry of the unit cell. To get a fast convergence, it 
is imperative that the sides of the triclinic unit cell are labelled such that 7 m i n is as large as 
possible. 
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Now, we consider the sum G20 defined in Eq. (J3.22|) . Using the relation, 

5 = 

which is derived in Appendix B, we can write 



4tT% 



G 20 (r) 



h 



E 

ni 



exp 


— 2irnit3 


6 




ni| 


6 





exp (27rm 1 Xi) — 



2tt U 



I, 



cosh ( 2nt 3 5 ) — cos {2ttxi) 



In (2) 



where 



(3.26) 



(3.27) 



Xi — — C 31 t 3 + ti, X2 — —032^3 + t 



2- 



(3.28) 



and we have used the identity from Eq. (j2.13|) . Thus, we have been able to obtain G2o( r ) 
analytically. Now, we transform G' 2 (r). The sum over n\ in Eq. (j3.21|) , 



S(n 2 ,n,r 3 ) = 

ni=— 00 

can be transformed using an identity, 



exp (27rmiXi) exp(-27T7 niin2 t 3 ) 



7m 



"2 



(3.29) 



exp -/3\/ a 2 + (q + nS)' 



exp [«p (q + nS)] 



(3.30) 



\J a 2 + (q + n5) 2 

=i|x>° (V^+H-^j exp ( 2?r v g ) ' 

which can be derived with a simple application of Jacobi Poisson theorem^ to the integral 



/ \ 1 f+°° exp (-b^ a 2 + y 2 

K (aVb 2 + x 2 j=- / dy / n = ^ exp (ixy) . 



a 2 + 



(3.31) 



Identifying 



X\ 

(3 = 2nt 3 , p = 2tc—, q = n 2 a, a 
5 



n 2 b 



and 5 = 5, 



[3.32) 
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one obtains 



S (n 2 , n, r 3 ) = exp ( -2m^n 2 a V" 



exp —2nt; : 



n 2 b 



+ (?2 2 a + ni<5 



n 2 6 



n 2 a + ni<5 



x exp 



exp 



( — 2Tti^n 2 a ) x K I 2ir 



.ni 



x exp 27ci—n 2 a . 
\ 5 J 

Substituting the value of S (n 2 ,ri,r 3 ) in Eq. (j3.21|) we obtain 



n\ — x\ 



(3.33) 



exp 



"2 



X)^o 2tt 



n 2 6 







7/. 




/m - 




V 5 



exp ( 27ri^n 2 a j . 



(3.34) 



Combining Eqs. (j3.20j) . f)3.22j) and ()3.34|) we get one of the main results of this paper, 



exp 



"2 



. , x x a 
2mn 2 x 2 — 

5 



(3.35) 



x 



J2 K ° \ 2rr 



n 2 b 



4 



ni 
1 



Hi — X\ 

6 



exp I 2ni^n 2 a 
5 



— — In | cosh ( 2nt 3 5 ) — cos {2-kx\ ) 



In (2) 



Result in Eq. 1)3.35)1 represents the sum for slab geometry and generalizes the results of Arnold 
et alM. Similar expressions were obtained by Liem et alr^. Substitution of G 2 from Eq. 1)3.35)1 
in Eq. (J3.14j) gives us an alternative form of G. We note that the problem of convergence 
with Eq. ()3.35)) still persists if the charges are close together. The slowest converging term 
in Eq. (j3.35j) goes as K (2n \n 2 \ p) and it still does not converge fast enough when, 



P = b\ tl+ ^ 



(3.36) 
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is small. The problem of convergence lies with only those terms corresponding to ri\ = 0. 
So, we separate out these terms 



where 



and 



G 2 — G 2 + G\d, 



(3.37) 



G'i = y exp 



lln 



2717712 X 2 - 



X 



n 2 b 



ni 



Xid\ 




8 )_ 




(U\- 


X\ 


\ 5 



(3.38) 



exp ( 2ni-^-n 2 a 



Gid = t y2 ex p 

h — 

"2 



27iin 2 ( x 2 



n 2 =l 



COS 



2im 2 x 2 - 



x±a 
8 

x{a 
8 



Kn 27T 



n 2 b 



(3.39) 



K I 27m 2 b\lt 2 3 + 



The term G 2 ' does not have any convergence problem for small separation between the two 
charges. We need to apply a final transformation to the sum in Eq. ()3.39|) . We start with 
the identity^, 



/ (p, x) = 4 K (2imip) cos (2 

m=l 

-{, + !»(§)} 1 



7rmx 



(3.40) 



iV-1 



p 2 + x 2 
+ 



Tit =1 



\J p 2 + {n\ + x) 2 \j P 2 + ( n i - X Y 
2 7 - {^(iV + ar) + - x)} 



+ £ 



-1/2 
/ 



p M (C (2Z + 1, iV + x) + ( (21 + 1,N — x)) 



where if) and ( stand for digamma and Hurwitz Zeta function respectively. iV > 1 is the 
smallest integer chosen such that it satisfies the condition N > p + x. However for better 
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convergence it is desirable that one chooses N such that N > p + 1. Now, identifying p from 
Eq.flirHl) and 

(3.41) 



x 



X 2 



6 



and realizing that (see Appendix B) 



o 9 r? + rl + r2 + 2rir 2 cos a + 2r 2 r 3 cos 3 + 2r 3 ri cos 7 
p + x = 2 



/ 2 
'2 



(3.42) 



one obtains an expression for G, which converges exponentially fast even for small xf. 



5 ' 



G(r) = Y~ 22 ex P [27ri(niti + n 2 t 2 )] B(ni, n 2 r 3 , £) 



711,712 



+ 



J^exp 



27rm 2 x 2 - 



5 



x K \2tt 
1 



n 2 b 



■tl + ' ni - Xl 



In 



+ 



7T 



Z 2 3 



cosh ^27rt 3 (5j — cos (27ra;i 

(i + «i) + %(0- 



exp 2ixi—n 2 a 
6 

In (2) 2 7 



2, /p\ t/j(N + x) + ip(N-x) 
In ( 

2 



/9 



1 OO 
71=1 

4e 



■!/2 I ,„ 



p 2n [C (2n + 1, N + x) + C (2n + 1, iV - x)] 



711=1 



\J p 2 + (ni + x) 2 \j P 2 + ( n i - x)' 



+ 



(3.43) 



(rj + r\ + rf + 2r^ 2 cos a + 2r 2 r 3 cos 8 + 2r 3 r x cos 7) 1 ^ 2 

Even though Eq. ()3.43|) gives a very good convergence for smaller values of < e = 10~ 3 , 
it is not defined when t 3 = and x\ = 0. The problem lies in the logarithmic terms, which 
can be combined together such that the opposing logarithmic divergences cancel each other 
as shown below. For small separations, it can be easily shown 7 that 



In [cosh y — cos x] = In 



y 2 + x 2 



+ In { 1 + § (y 2 - x 2 ) + I (y 4 - x 2 y 2 + x 4 ) (3.44) 



+| (V 4 + ^) (y 2 -x 2 ) + 0[x 8 ,y*]y 
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Thus all of the logarithmic terms in Eq. (j3.43|) can be combined together as 
1 



1 2 



cosh ( 2irts5 ) — cos (2irx 



(3.45) 



2, /4tt5\ , ( 2! 
'fe IT/ + { 4! 



US 



.} I 



2! 
+ 6! 



US) -x\ (us) 2 + xi 



2! 

+ +777 



US) +x\ 



us 



x. 



+ 



x*, [us 



The RHS of Eq. ([3.45|) remains regular even when x\ and U both tend to zero. The 
self-energy of the system can be easily obtained now as, 



G s 3 e »= lim I G( ri ,r 2 ,r 3 ) 

(ri,r2,r 3 )— »(0,0,0) 



(rf + r\ + rf + 2rxr 2 cos a + 2r 2 r 3 cos (3 + 2r 3 ri cos 7) 

(3.46) 



1/2 



^ ^ s(m, 712,0,0 



ni,ri2 
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x exp ( 2mn\n2- 



8 



\ 5 \ tt 2 f 4ttS\ 27 
~kS ~T 2 MX/ 



We have thus obtained complete expressions for G and the self-energy. 



IV. RESULTS AND CONCLUSIONS 



We have obtained complete expressions for the logarithmic potential in 2D and Coulomb 
potential in 3D, including the self-energies. The results were derived for most general cases, 
that is a rhombic cell in 2D and a triclinic cell in 3D. To my knowledge, this is the first time 
a practical method has been developed in 3D, which is different from the Ewald method, 
and yet may be applied to a triclinic unit cell to obtain periodic Coulomb sums. Even 
though the formulas developed here look complicated, their implementation on a computer 
will be marginally difficult from the case of orthorhombic unit cell. The formulas derived 
here converge extremely fast and require only a few dozen terms at worst to obtain results 
to a very high accuracy as opposed to the Ewald method, which may require close to 200 to 
300 terms for the same calculations. In the process, we have simplified and solved a problem 
mentioned by Crandallifi, that of finding Coulomb potential in close vicinity of a particle 
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under the PBC An important implication of the formulas derived here is that most part of 
the interaction can be calculated linearly in the number of charges present in the system. 
For more details on how this can be achieved, we refer the reader to Sperb who discusses 
this in the context of an orthorhombic cell. The results obtained in this paper reduce to the 
results of a recent paper— when all angles pertaining to the unit cell are set to 7r/2. 

The results for 3D triclinic case may be obtained by directly generalizing Lekner's work. 
This work by the author will be presented elsewhere. We also note that the logarithmic 
sum in 2D for a rhombic cell may be obtained in a closed form. This will be the subject 
of another paper. Also, here we would like to point out a connection between the results 
of slab geometry and that of 3D triclinic cell. As it has been shown here, the 3D Green 
function can be broken in two parts. The first part corresponds to the slab geometry Green 
function and the second part takes into account the rest of the layers. Thus following Refll6 



one can make use of this relation to obtain potential energies for the slab geometry cases by 
employing the result for the triclinic cell. 

A naive application of most methods gives a scaling which goes as N 2 , where N is the 
number of charges in the unit cell. However, the Ewald method can be optimized^ 2 - to give 
a scaling of [N ln(iV)] 3 / 2 . Strebel et aZ.— gave an approximate method, which they call the 
MMM method, which is based on the formulas developed by Sperb in his earlier work 2 -. With 
the help of MMM one can achieve a Nln(N) scaling. There is another approximate method 
in use to achieve a faster scaling. This method is known as PPPM. In Refllll however, it 
was shown that for iV > 2 10 and a relative tolerance of 10~ 4 the MMM is the best method 
available. As the formulas derived here are a generalization of Sperb's work, it may now be 
possible, by using the results presented in this work, to employ the MMM method to achieve 
a scaling of iV ln(iV) even for a triclinic cell. Similarly for the logarithmic interaction in 2D, 
it should be possible to achieve N\n(N) scaling. 

In short, we believe the method developed here is an alternative to the Ewald method. 
The formulas developed here generalize and simplify Sperb's 2 work. From the results of 
summation formula derived for a triclinic cell, it is easy to obtain results for 2D + h slab 
geometry and vice versa. For the slab geometry case expressions obtained here generalize 
the work of Arnold et al lA . and give an alternate derivation of the results obtained by Liem 
et a/.—. The formulas derived in this work can be easily employed to calculate the Madelung 
potential for any periodic crystal in 3D, where a triclinic cell repeats itself to infinity under 
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(N+1/2,-y) 



FIG. 2: Contour of integration, Cn- Both N and y tend to infinity. 



the PBC. 
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APPENDIX A: COMPLEX SUM 

The usual way to sum over series of type 

oo 

S= £ f(n) 

n=— oo 

is to consider the integral 



(Al) 



I = j f (z) 7T cot (nz) dz. (A2) 

It is required that the function / (z) satisfies the condition that integral I becomes zero 
when the contour of integration is chosen to be Cm as shown in Fig. 2. 
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The poles of 7rcot (ttz) fall at z = n where n = 0, ±1, ±2, .... Then by residue theorem 
we have 



2_j f ( n ) = ~ sum °f residues of / (z) n cot (ttz) at the poles of / (z) . (A3) 

n=— oo 

Here, in particular, we consider the function 

„ . , exp (ina) n . , A , 

f n = : L.a « < 27r ' ( A4 ) 

where x, f3, 7 are real numbers and x > 0. The results obtained here are more general in 
nature and may be applied for other forms of f(n). It can be easily verified that the function 
given above does not satisfy the condition that integral / go to zero for contour CV. A trick 
which is usually not found in books may help solve the problem. Instead of considering the 
integral / in Eq. (jA2|) . we consider the following integral 

I' = £f(z) (-l) z n csc (nz) dz, (A5) 

where have in mind that exp (— wr) = —1. Residues of / (z) (— l) z 7r csc (ttz) at z — n, 
n = 0,±1,±2, is 

lim / (z) (z - n) (-1)*tt csc (ttz) = f (n) . (A6) 

z— *n 

Thus if the integral /' goes to zero for the contour CV then we obtain 

/ (n) = —sum of residues of / (z) (— l) z n csc (ttz) at the poles of / (z) . (A7) 



Function / given in Eq. (jA4j) does satisfy the condition that I' = when the integration is 
evaluated for the contour Cat. To show this, we concentrate on the the following function 

. . exp (izat) exp (— inz) , . . 

9 (x, y, a) = / H V '-. A8 

sin (ttz) 

We note that in terms of g (x, y, a), f(z) can be written as 

m = , ; ■ (A9) 

(n — (3) + Y 



With the substitution of z = x + iy in Eq. (jA8|) . we obtain 



2 exp (—yea) 
[exp (— Any) — 2 cos (2irx) exp (— 27vy) + 1] J 



(x, y,a)\ = — — - — — — — ~ /2 . (A10) 
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0.50 



FIG. 3: The value of function <? max remains between 0.5 and 1.0 as a is varied from to 2tt. 



We note that 



lim \g(x,y,a) \ = 0. 

|]/|-+0O 



(All) 



The condition in Eq. (jAllj) ensures that I' goes to zero on those portions of the contour 
which lie parallel to the x axis. To consider the portions of contour parallel to the y axis, 
we substitute x = N + 1/2. One obtains 

2 exp (—yat) 



\g(N + l/2,y,a)\-. , 

1 + exp (— ziry) 

which implies that the maximum value of the function \g (N + 1/2, y, a) | occurs for 

1 , (In -a 

y = 7r ln \ 



(A12) 



and is given by 



2tt — a 



a 



2tt 



exp 



a . f 2n — a 

In 

2tt V a 



(A13) 



(A14) 



A plot of g max (a) is shown in Fig.3. 

It is clear that the value of g max (ct) remains between 0.5 and 1.0, and this ensures |/(z)z| 
goes to zero on the contours parallel to the y axis. Thus it is clear that I' goes to zero when 
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evaluated for the contour Cn and hence by the application of formula in Eq. (jA7|) we obtain 

E exp (ma) exp \(/3 + ry) (a + in)] r . 
, 0,2 ' 2 = 9 ' csc W + «7)] (A15) 

exp[(/3-z7) (a + ivr)] 

2^ £ 7T CSC [TT{p — t'-fJl 

ir exp [i/3 (a — 2ir)] sinh (7a) + exp (i(3a) sinh [7 (2tt — a)] 
7 cosh (2717) — cos (27r/?) 

APPENDIX B: TRICLINIC CELL 

Let us consider the most general type of triclinic cell shown in the Fig. 4. The cell is 
characterized by sides li, l 2 and Z 3 and angles a, (3 and 7. We can choose the unit vectors 
along the directions of the triclinic cell as 



d = (1,0,0), 

e 2 = (cos a, sin a, 0) 

cos f3 — cos a cos 7 



(Bl) 



e 3 = cos 7, 



sin a 



sin 2 7 



cos (5 — cos a cos 7 



sin a 



We can now get reciprocal vectors using equation. Now we can calculate c^- and we get the 
following results using the package Mathematica 



_ Z3 sin 2 (3 sin 2 7 _ 

Cn ~~ 72^~2 — ' C 22 — 72 • 2 I ' C 33 ~ l^ij 



Cl2 



<f sm a <2 sm a 

'3 



i| / — cos a + cos /3 cos 7 



Z1Z2 V sin a 



l 3 ( — cos (3 + cos a cos 7 

C23 = T 



h \ sin a 

I3 { — cos 7 + cos a cos /? 

Cl3 = T 



h V sin 2 a 

We can then obtain p and x as follows 

(r 2 cos 2 a + rf cos 2 /? + 2rir 3 x [cos 7 — cos a cos f3\ ) ^ 2 
P = j ( B3 ) 

and 
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FIG. 4: A triclinic cell explaing different labels for sides and angles. 



x 



r 2 + r\ cos a + r 3 cos (3 
To ' 



Finally we obtain 



2 2 r\ + r\ + rf + 2r!r 2 cos a + 2r 2 r 3 cos/3 + 2r 3 r x cos 7 
p ~\~ x — 



''2 



Using the relations given above, it can be shown on Mathematica that 



(B4) 



(B5) 



Vbl 



4vr 2 / 



2- 



(B6) 



where 63, V and 5 are defined in Eqs. (j3.5|) . f)3.2jl and f)3.24j) . 
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